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Abstract. The denaturation dynamics of a long double-stranded DNA is studied 
by means of a model of the Poland-Scheraga type. We note that the linking of the 
two strands is a locally conserved quantity, hence we introduce local updates that 
respect this symmetry. Linking dissipation via untwist is allowed only at the two ends 
of the double strand. The result is a slow denaturation characterized by two time 
scales that depend on the chain length L. In a regime up to a first characteristic 
time Ti ~ L^'^^ the chain embodies an increasing number of small bubbles. Then, 
in a second regime, bubbles coalesce and form entropic barriers that effectively trap 
residual double-stranded segments within the chain, slowing down the relaxation to 
fully molten configurations, which takes place at T2 ~ ■ This scenario is different 
from the picture in which the helical constraints are neglected. 



PACS numbers: 87.14.gk, 87.15.H-, 61.25.hp, 36.20.-r 



Multiple timescales in a model for DNA denaturation dynamics 



2 



The Watson-Crick helix is the typical form of DNA in the cell pQ. In the laboratory, 
upon heating DNA molecules in solution, one obtains an helix-coil transition called 
denaturation. Since decades, DNA denaturation has attracted the attention of scientists 
because it can help to understand important biological processes: for instance, the 
genetic code can be accessed during transcription and replication by an opening of 
bubbles [1]. It is experimentally known that the fraction of molten DNA increases for 
increasing temperature [2j. The first theoretical description of denaturation that can 
account for this phenomenon was a simple model by Poland and Scheraga (PS) [3]. 
A model of this kind is now behind software like Meltsim [4j, predicting sequence- 
dependent melting curves, which can then be compared with experiments. Another 
simple model for DNA is due to Peyrard and Bishop (PB) [5], [6]. Also some more 
detailed yet mesoscopic models have been recently proposed [3 [8], allowing for the 
numerical study of features that cannot be simulated in all-atoms molecular dynamics. 

PS models rarely take into account the topological state of a macromolecule. All 
polymers that form closed rings have some conserved topological features, as long as 
the chains cannot break and cross each other. For instance, in circular DNA, such 
as the genomes of some bacteria, the number of times that the two strands twist 
around each other (linking number) cannot change. If this constraint is included in 
a PS model, the thermodynamic denaturation transition is essentially suppressed [9], 
unless supercoiling effects are considered [10]. The topology is not fundamental for 
the equilibrium properties of linear polymers, but the fact that chains cannot cross each 
other is clearly relevant for their dynamics. For instance, in electrophoresis the dynamics 
of DNA under an electric field depends on its continuous entangling with the polymers 
of a gel [H]. This feature is taken into account in models of electrophoresis [H], like 
the Rubinstein-Duke "repton" model [T2| [T3| fll] . 

While the equilibrium thermodynamics of DNA has been largely investigated, 
the dynamics of this macromolecule is of interest as well. The dynamics of thermal 
denaturation has been studied by means of the PS and similar models [151 IIHl [13 [T8| - 
by using the PB model [6], and in more detailed models [TIE]- For short DNA segments, 
the rates of opening found in experiments [19] can be estimated by PS models with 
stochastic dynamics [17]. However, being a coarse grained description, the PS model is 
particularly useful to study long chains. 

The aim of this paper is the study of the denaturation dynamics of a long DNA in a 
PS model with a stochastic evolution that takes into account the helical structure of the 
double strand. This was not explicitly included in the standard formulation of the PS 
model [iSl HE] . We show that a dynamics with local preservation of the linking between 
the chains yields a new scenario, with two time scales. A first regime is dominated 
by denatured bubbles diffusing into the chain from its ends (where the double chain 
can freely untwist). This regime ends after a time lapse that scales approximately 
as the square of the chain length, where the number of bubbles reaches a maximum. 
Then bubbles start to coalesce, with the ones near the chain ends that form entropic 
barriers trapping the helical domains still in excess. The trapping into these metastable 
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states further slows down the denaturation process in the model, leading to a thermally 
equilibrated denaturation only after a second time scale, which grows as the cube of the 
chain length. 

In the PS model, each strand is represented by a chain of length L, where a site i 
(1 < i < L) stands for a local portion of DNA. We associate each site to a segment of 
10 base pairs, which thus represents a complete helical turn when paired 0. The state 
of the chain is stored in a Boolean array (Xj, where o"i = 1 if the two segments at index i 
are paired and o"j = otherwise [see Fig. [T]^a)]. DNA conformation are thus represented 
by an alternation of segments of paired bases (sequences of I's) and of open bubbles 
(sequences of O's). For every a,, = 1 there is a binding energy e = — 1. At a temperature 
T this corresponds to a Boltzmann factor q = e"*^/^. Thus, a sequence of m paired 
bases brings a contribution to the global weight W of the configuration. A bubble 
formed by two complementary strands of length i instead has an entropic contribution 
accounting for all the possible conformations of a walk of length 2t. if s is the entropy 
per step of a walk, the constraint to form a loop yields a weight Bs'^^i~'^, where B is 
a constant factor. The exponent c can be deduced from self-avoiding walks statistics: 
with the excluded volume between the chains fully taken into account [20] it is c ~ 2.1. 
Hence, the weight of a whole configuration is 

W = {q"''){Bs^^H^^) . . . {q"'^){Bs'^^"i-''){q'''^+') (1) 

where u is the number of bubbles. We also set ctq = (Jl+i = 1, namely each end of a 
single strand is joined to the corresponding end of the other strand (no Y-fork is formed). 
Thus, at high T the equilibrium configuration is a ring of length 2(L + 1). At low T, 
on the other hand, typically one finds long double-stranded parts separated by small 
bubbles. According to this description, the properties of the model at thermodynamic 
equilibrium can be derived analytically [3l [20] . 

The simplest dynamical rules that can be assigned to the PS model involve moves 
where locally one cTj changes, 

ai = l < — > ai = 0. (2) 

A Metropolis criterion can then be used to choose whether to accept the move. This 
kind of update resembles the dynamics of adsorption of a polymer onto a wall, see 
Fig. [Tt^a)-(b). However, for 1 <^ i <C L, we note that the (dis) appearance of a "1" would 
imply either a temporary breaking of one of the chains to (un)twist the two strands 
there (as it happens e.g. with topoisomerase enzymes [21]) or a global rotation of 27r 
of the whole part < i of the chain with respect to the whole part > i. The latter case 
is not in agreement with the idea of small time step that is intrinsic in ([2]). Since the 
update ([2]) neglects the helix (and the consequent link) of the DNA strands, one needs 
another dynamics that respects the local topology of dsDNA. 

I Note that this choice of the coarse graining is not fundamental, as we are just interested in the scahng 
properties of long chains and not in the microscopic details of the dynamics. 
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In order to preserve locally the linking number, we adopt a different basic move: 
one piclfs a boundary {i\i + 1) at random and swaps the relative variables, 

(Ji = X, (Ji+i = y — > ai = y, (Xi+i = x , (3) 

where x and y can be or 1 (if they are equal, the move is trivially the identity). This 
exchanges the amount of linking of the chains at position i with that at position i + 1, see 
the sketch in Fig. [D^c). In practice, the linking plays here the role of an order parameter 
that is locally conserved. In order to have an evolution of the global linking number, 
however, one has to consider also some dissipation via untwisting of the double helix. In 
fact, untwisting with an update like ([2]) should be valid close to the ends of the dsDNA 
macromolecule, if they are free to rotate. Two special updates are then introduced by 
adding to the list of possible boundaries the (0|1) and the (L|L + 1) ones: 

0"! — > 1 — (Ji if (0|1) is chosen (4) 
^ 1 — (Tl if (L|L + 1) is chosen (5) 

These two moves can thus change the energy E = cij of the chain, allowing 
equilibration at every temperature. It is important to note that the move ([3]) can 
lead to the nucleation of new bubbles along the whole chain, from the boundaries of 
already present ones (e.g. a = . . . 00001111 cr = .. . 00010111 . . .). The dynamics 

obeys detailed balance, hence clearly also the reverse process of bubble coalescence 
can take place. Therefore, move ([3]) seems the best approximation for the purpose 
of describing the local conservation of the linking between two complementary DNA 
chains. By definition this model cannot deal with the formation of twisted bubbles, 
for example by means of a breaking of base-pair bonds without untwisting the chain. 
These configurations, however, are entropically unfavored compared to the ones in which 
the twist is concentrated on paired segments and open bubbles are expanded. We 
thus expect that the approximation of neglecting the formation of twisted bubbles is 
appropriate in a simple model. 

A time step consists in a sequence of L realizations of a basic move, each one with i 
picked at random, uniformly along the chain. If the dynamics ([3])-([5]) is used, i G [0, L], 
while i G [1,-^^] for update ([2]). Then, according to the Metropolis criterion each move 
is accepted with probability p = min{l, PVnew/W^ow}, where PFnew is the weight of the 
proposed configuration and Wo\d is the weight of the present configuration. 

The protocol on which we focus is a quench of a system equilibrated at low T to a 
regime at very high T. This is indeed the regime that allows us to appreciate more the 
effects of the new physical ingredients in this model. Equilibrated configurations to start 
the protocol are generated by setting g/s^ = 100 and by applying multiple ([2]) updates[§| 
[this because they equilibrate faster than ([3])]. Then, each protocol starts by switching 
instantaneously to q/s"^ = 0.01 at time t = 0. We set the parameter B = 1, postponing 
the systematic study of different cases to future works. For the bubble exponent we use 
the value c = 2.14 |22|. 



§ Similar results can be obtained by starting from the fully ordered state, ct^ = 1 for all i. 
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Figure 1. (Color online) (a) Sketch of a DNA configuration and of the relative array 
a. (b) Configuration obtained by updating (a) with ([2]) at site i. (c) Configuration 
obtained by updating (a) with ([3]) at the same site: in this case there is a local 
conservation of the linking number. 

The transient to the new equihbrium is monitored by studying the scahng properties 
of two quantities, the number of denatured pairs 5 = ~ ^i) number 

of bubbles v. The former is the quantity normally inferred from UV- absorption 
experiments [2| while v is useful for characterizing the state of the system. For 
convenience, data are binned in time intervals with constant size = log 1.05 in log-scale. 
Moreover, an average over at least 1000 trajectories is perfomed. Fig. [2] shows 5 and p 
vs time in log-log scale for L = 500, both for a dynamics involving only move ([2]) and 
for the dynamics ([3]) introduced in this paper. In the former fast denaturation 

takes place, at a time scale tq ~ 1 that does not scale with the system size L \1Q\ . 

The dynamics ©-([S]) generates a richer picture. Two characteristic time scales 
have been highlighted by vertical lines in Fig. [21 The process goes as follows: at high T, 
the rate di = 1 o"i = is much higher then the rate of the opposite transition. The 
same is true for site i = L. Thus, O's enter at the boundaries and diffuse toward the 
center of the chain. This first regime is characterized by an increase of 5 and u that is 
slower than ~ t^^^ [jj], see Fig. |2l A time scale Ti is characterized by the maximum of the 
number of bubbles and it marks the end of the first regime. At times t > ti one observes 
a second regime in which 6 continues to increase while u decreases, which implies that 

II In this regime we do not observe simple dynamical scaling for 6 and v. 




Figure 2. (Color online) Log-log plot of the number of open sites 6 (dense lines) and 
of the number of bubbles v (dashed lines) vs time, for L — 500. Data shown with thick 
lines are obtained with our update ©-([S]), while thin lines (red online) correspond to 
data obtained with update ([2]). The final state for both dynamics is the equilibrium 
at high T, with S ~ L and « 1. The dot-dashed line represents a scaling ^ t^^^. 
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Figure 3. (Color online) (a) Number of bubbles as a function of time, for chain lengths 
L = 100, 200, 300, 400, 500, 600, 800 and 1000 (from bottom to top), (b) Number of 
open sites vs t, same notation. 



bubbles coalesce. Finally, equilibrium is reached at a time T2, when 6 ^ L and u ^ 1. 

Both Ti and T2 increase with the system size, see Fig. |3l We find that both time 
scales are consistent with an algebraic dependence on L, i.e. 

Ti ~ , with zi ~ 2.15 ± 0.10 (6) 
T2 ~ , with Z2 ~ 3.0 ± 0.1 (7) 

The peak of the u plots is a particularly clear feature that helps to estimate the value 
of Zi. peaks are reached at ti ^ 0.04 x L"^-^^. Moreover, around ti we achieve a data 
collapse of the form u/L vs t/L^^: Fig. IH^a) shows that a critical density of bubbles 
v{ti)/ L is reached at ti. The values of the critical density of bubbles, as a function of 
L^^/^ and extrapolated for L 00, tend to 0.280(1). A similar collapse S/L"'^ vs t/L^'^ 
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Figure 4. (Color online) Curves rescaled to collapse at ti with Zi — 2.15, and at T2 
with Z2 ~ 3, see the notation of Fig. [3] and Eqs. (I6])-([7|. 

can be attained: the exponent yielding the best rescahng is a\ ~ 0.94, see Fig. |l](b). 

Data collapses can be done also to determine T2- At T2 by definition we have the 
full denaturation, i.e. 5{t2) ~ L. To estimate 2:2 ~ 3 we have in fact required that 
S{t)/L — >• 1 for t T2 (for all L's for which T2 could be reached by simulations), as 
shown in Fig. |l](d). A data collapse of z//L"2 vs t/L^, with 02 = 0-8 [Fig- 111(c)], confirms 
that Z2 — 3 is the exponent characterizing the time scale T2. 

The first regime is essentially a diffusion of random walkers, the CTj = entering 
from the boundary, and one should expect a temporal domain scaling as the square 
of the system size. Indeed, we estimate 2:1 = 2 for c = 0, a case in which the open 
sites are independent of each other and the weight of a configuration just depends on 
6. For c = 2.14 we instead estimate the small deviation zi ^ 2.15 from this classical 
result, probably due to the bubble weights. The case c = is also interesting because 
it does not display two different time scales but only one. It confirms that the bubble 
interaction and coalescence is the process leading to that second time scale. 

It is possible to predict the value of Z2 with a simple argument: in the regime 
between ti and T2 two bubbles at the chain ends trap the double-stranded parts inside 
the system, preventing a fast escape of (jj = 1 from the boundaries. Let us concentrate 
on one of the two ends, say i = 1, as shown in Fig. [5l where an exemplified escape of 
a double-stranded segment is shown. Taking the length £2 of the forming loop on the 
right as a reaction coordinate, the free-energy F = ~lnW has a profile like the one 
shown in Fig [5l It achieves a maximum at state (c), where £2 = ^i/2. The rate of 
escape from (a) to (e) is proportional to the barrier jump rate, which is proportional 
to the ratio of the weights in (c) and (a), [(£i/2)~'^]^/£^'^ ~ £]^^. Hence the time spent 
for this escape scales as {£iY- As this has to take place for all £1 up to the system size, 
T2 ~ ^i^^ii^iY ~ L'^^^, which would imply Z2 = c+ 1 ~ 3.14 if the whole dynamics was 
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Figure 5. (Color online) Sketch of the exemplified process of escape of double-stranded 
segments. Snapshots of some intermediate states and the corresponding free-energy 
profile (as a function of the length ^2 of the growing loop) are shown. 

as simple as the one described. Of course, multiple pathways intersect and a description 
only in terms of a single li might not be exhaustive. Nevertheless, this prediction is 
remarkably close to the value ^2 — 3 obtained by data collapse, suggesting that this is 
the main mechanism acting in the second regime. 

Besides time scales, we have also a time dependence of the number of open bases 
6 that is different from the one predicted by previous models, where one can observe 
a linear increase of 5 with time [IHl E] or 5 ~ t^/^ [T5], in conditions similar to the 
one discussed in this paper (initial state at low T, denaturation at high T, no external 
forces). In our model, on the other hand, we observe that 5 increases slower than the 
square root of time. 

The fact that polymer chains cannot cross each other is at the basis of our version 
of the PS model, but of course it is included in many other models, like the model of 
polymer diffusing into a gel by de Gennes [23], in which he found that the diffusion 
constant scales as Furthermore, in simulations of polymers in dense melts [21] 

one observes autocorrelation times scaling as L^. These long time scales derive from 
the reptation dynamics of the polymers, which have to diffuse into tubes formed by the 
melt. We argue that scenarios like that one are similar to the phenomenon predicted by 
our DNA model, because the two twisted strands constrain the stochastic movements 
of each other in space during melting. 
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